phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p);

dimensionedScalar myNu
(
    "1",
    dimTime * dimensionSet(0, 2, -2, 0, 0) / p.dimensions(),
    1
);

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
    fvScalarMatrix pEqn
    (
        fvm::laplacian( myNu, p ) == fvc::div(phi)
    );

    pEqn.setReference(pRefCell, pRefValue);
    pEqn.solve();

    if (nonOrth == nNonOrthCorr)
    {
        phi -= pEqn.flux();
    }
}

Info<< "continuity error = "
    << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
    << endl;

U = fvc::reconstruct(phi);
U.correctBoundaryConditions();

Info<< "Interpolated U error = "
    << (sqrt(sum(sqr((fvc::interpolate(U) & mesh.Sf()) - phi)))
      /sum(mesh.magSf())).value()
    << endl;
